home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YCOMDE.C < prev    next >
C/C++ Source or Header  |  1988-11-28  |  6KB  |  226 lines
  1. /*******************************  YCOMDE.C   *********************************/
  2. /***************** Differential Equation COMMANDS AND MENU *******************/
  3. /********************* (C) 1986,7,8 by JAMES A. YORKE ************************/
  4.  
  5.  
  6. #include "yinclud.h"
  7.  
  8. /* includes invert() */
  9.  
  10. static int(*savePlotPointer) () = NULL;
  11.  
  12. int     DECommands(CodeName)
  13. register char    *CodeName;
  14. {
  15.     int     st;
  16.     int     connectp();    /* needed for pointers */
  17.     int     plot(), euler(), rungekutta(), iterate_map();
  18.                 /* needed for pointers */
  19.  
  20.     TEST2("conn", "connect") {
  21.         plotFlag = YES;
  22.         if(SCREEN)
  23.             PRINT "Now connects dots with lines.    \n");
  24.         savePlotPointer = plotPointer;
  25.         plotPointer = connectp;
  26.         return(1);
  27.     }
  28.  
  29.     TEST2("dis", "disconnect") {
  30.         plotFlag = 0;
  31.         if(SCREEN)
  32.             PRINT "Now plots disconnected dots.    \n");
  33.         if(savePlotPointer != 0L)
  34.             plotPointer = savePlotPointer;
  35.         savePlotPointer = 0L;
  36.         return(1);
  37.     }
  38.     TEST("e") {
  39.         if(step != -9999.) {
  40.             PRINT
  41.                 " E: computes single step error every 5 steps; (toggle); currently");
  42.             ODEcheck = 1 - ODEcheck;
  43.             toggle(ODEcheck);/* prints " ON\n" or " OFF\n" */
  44.             return(1);
  45.         }
  46.     }
  47.  
  48.     TEST("euler") {
  49.         IterIteratee = euler;
  50.         return(1);
  51.     }
  52.     TEST2("v", "invert") {
  53.         invert();
  54.         return(1);
  55.     }
  56.  
  57.     TEST("runge") {
  58.         IterIteratee = rungekutta;
  59.         return(1);
  60.     }
  61.  
  62.     TEST("step") {
  63.         step = Entervalue(step, CHECKSET);
  64.         sixth_step = step / 6.;     /* for rungekutta */
  65.         halfstep = step / 2;
  66.         return(1);
  67.     }
  68.  
  69.     TEST("spc") {
  70.         if(steps_per_cycle != -9999) {
  71.             if(SCREEN) {
  72.                 PRINT
  73.                     "\n");
  74.                 PRINT
  75.                     "This is the number of differential equation steps in one period of the \n");
  76.                 PRINT
  77.                     "of the forcing period.\n");
  78.                 PRINT
  79.                     "     ******IT MUST NOT BE SMALL******* \n");
  80.                 PRINT
  81.                     "The larger the value, the more accurate the solution.\n\n");
  82.                 PRINT
  83.                     "See command \"IPP\" (in menu \"P\") to change the plotting frequency.\n\n");
  84.             }
  85.             st = steps_per_cycle;
  86.             steps_per_cycle = Entervalue(steps_per_cycle, CHECKSET);
  87.             step = (twopi / C4) / steps_per_cycle;
  88.             sixth_step = step / 6.;/* for rungekutta */
  89.             halfstep = step / 2;
  90.             if(its_per_plot == st)
  91.                 its_per_plot = steps_per_cycle;
  92.             if(steps_per_cycle < 10 &&
  93.                     steps_per_cycle > -10)
  94.                 PRINT
  95.                     "WARNING: THIS IS PROBABLY TOO SMALL. YOU NEED LOTS OF STEPS PER CYCLE. \n\n");
  96.         }
  97.         else {
  98.             erase_line();
  99.             PRINT "Improper category for this map \n");
  100.         }
  101.         return(1);
  102.     }
  103.     return(0);
  104. }
  105.  
  106.  
  107. DEMenu() {
  108.     if(level == SETPARAM)
  109.         scr_clr();    /* in desmets pcio.a */
  110.     scr_rowcol(0, 2);
  111.  
  112.     PRINT
  113.         "                       DIFFERENTIAL EQUATION MENU                      \n\n");
  114.  
  115.     if(steps_per_cycle == -9999.&& step == -9999.)
  116.         PRINT
  117.             " The current process is not a differential equation\n\n\n");
  118.  
  119.     PRINT
  120.         " CONN: CONNects consecutive dots as plotted with a line segment       \n");
  121.     PRINT
  122.         " DIS:  DISconnect turns off CONNECT; plots dots, not line segments \n\n\n");
  123.     if(step != -9999.) {
  124.         if(steps_per_cycle == -9999.)
  125.             PRINT
  126.                 " STEP: Step size for differential equation = %g             \n", step);
  127.         if(steps_per_cycle != -9999)
  128.             PRINT
  129.                 " SPC: Steps_Per_Cycle(dif eq) %lf                      \n", steps_per_cycle);
  130.         PRINT
  131.             " EULER: Euler first order differential equation solver, fixed step size  \n");
  132.         PRINT
  133.             " RUNGE: Runge-Kutta differential equation solver, 4th order, fixed step  \n");
  134.         PRINT
  135.             " E: computes single step error every 5 steps; (toggle); currently");
  136.         toggle(ODEcheck);/* prints " ON\n" or " OFF\n" */
  137.  
  138.         if(max_error > 0)
  139.             PRINT "maximum dif equ error = %11.8le \n", max_error);
  140.         PRINT "\n");
  141.     }
  142.     PRINT
  143.         " I: inverts differential equations and the Henon map\n");
  144.  
  145.  
  146.     return;
  147. }
  148.  
  149.  
  150.  
  151. int     invert() {        /* 'v' inverts henon: map        x ->
  152.                    rho - x*x + C1*y inverse     -y/C1 ->
  153.                    rho/(C1**2) -(-y/C1)**2 +(-x/C1)/C1    */
  154.     int     i;
  155.  
  156.     if(strcmp("h", MapName) == 0
  157.             || strcmp("H", MapName) == 0) {/* if map == henon */
  158.         if(C1 == 0) {
  159.             PRINT "YOU CANNOT INVERT THIS BECAUSE C1 = 0   \n");
  160.             return(NO);
  161.         }
  162.         if(level >= PROCESS) {
  163.             turnoff(BIGCROSS);
  164.             turnoff(SMALLCROSS);
  165.         }
  166.         X_lower = -X_lower / C1;
  167.         Y_lower = -Y_lower / C1;
  168.         X_upper = -X_upper / C1;
  169.         Y_upper = -Y_upper / C1;
  170.  
  171.         flipHenon(y);
  172.         flipHenon(ya);
  173.         flipHenon(yb);
  174.         flipHenon(yc);
  175.         flipHenon(yd);
  176.         flipHenon(ye);
  177.         for(i = 0; i < 8; i++)/* y1,... y8 */
  178.             flipHenon(y + eqn0 + i * eqn1);
  179.         rho = rho / (C1 * C1);
  180.         C1 = 1 / C1;
  181.  
  182.         X_coord = 1 - X_coord;
  183.         Y_coord = 1 - Y_coord;
  184.         setXYPlot();
  185.         ScreenConstants();/* defines ax,bx where x = 0,00,1,2 */
  186.  
  187.         PRINT
  188.             "\n\n Map inverted; rho = %g; Jacobian C1 = %g\n"
  189.             ,rho, C1);
  190.         PRINT
  191.             "X_coord=%d;  Y_coord=%d                 \n", X_coord, Y_coord);
  192.         return(YES);
  193.     }
  194.     else
  195.         if(step != -9999.) {/* if we have a differential equation */
  196.             step = -step;
  197.             halfstep = -halfstep;
  198.             sixth_step = -sixth_step;
  199.             if(step < 0.)
  200.                 PRINT
  201.                     "INVERSE PROCESS: time step = %g                \n", step);
  202.             else
  203.                 PRINT
  204.                     "NOW BACK WITH DIRECT PROCESS: time step = %g                \n"
  205.                     ,step);
  206.             return(YES);
  207.         }
  208.     PRINT
  209.         "THIS ROUTINE(invert) DOES NOT HAVE THE INVERSE AVAILABLE. ABORT COMMAND\n\n"
  210.         );
  211.     return(NO);
  212. }
  213.  
  214.  
  215. flipHenon(yy)
  216. double *yy;
  217. {
  218.     double  temp;
  219.  
  220.     if(*yy == -9999.)
  221.         return;
  222.     temp = -yy[0] / C1;
  223.     yy[0] = -yy[1] / C1;
  224.     yy[1] = temp;
  225. }
  226.